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Oscillatory systems with time-delayed pulsatile feedback appear in various applied 
and theoretical research areas, and received a growing interest in recent years. For 
such systems, we report a remarkable scenario of destabilization of a periodic regular 
spiking regime. At the bifurcation point numerous regimes with non-equal interspike 
intervals emerge. We show that the number of the emerging, so-called “jittering” 
regimes grows exponentially with the delay valne. Although this appears as highly 
degenerate from a dynamical systems viewpoint, the “multi-jitter” bifurcation occurs 
robustly in a large class of systems. We observe it not only in a paradigmatic phase- 
reduced model, but also in a simulated Hodgkin-Huxley neuron model and in an 
experiment with an electronic circuit. 
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Interaction via pulse-like signals is important in neuron populations ra, biological mu, 
optical and optoelectronic systems [U]. Often, time delays are inevitable in such systems as 
a consequence of the hnite speed of pulse propagation [7] . In this letter we demonstrate that 
the pulsatile and delayed nature of interactions may lead to novel and unusual phenomena in 
a large class of systems. In particular, we explore oscillatory systems with pulsatile delayed 
feedback which exhibit periodic regular spiking (RS). We show that this RS regime may 
destabilize via a scenario in which a variety of higher-periodical regimes with non-equal 
interspike intervals (ISIs) emerge simultaneously. The number of the emergent, so-called 
“jittering” regimes grows exponentially as the delay increases. Therefore we adopt the term 
“multi-jitter” bifurcation. 

Usually, the simultaneous emergence of many different regimes is a sign of degeneracy 
and it is expected to occur generically only when additional symmetries are present [21E]. 
However, for the class of systems treated here no such symmetry is apparent. Nevertheless, 
the phenomenon can be reliably observed when just a single parameter, for example the 
delay, is varied. This means that the observed bifurcation has codimension one [9]. In 
addition to the theoretical analysis of a simple paradigmatic model, we provide numerical 
evidence for the occurrence of the multi-jitter bifurcation in a realistic neuronal model, as 
well as an experimental conhrmation in an electronic circuit. 

As a universal and simplest oscillatory spiking model in the absence of the feedback, we 
consider the phase oscillator dtp/dt = cn, where (p G M ( mod 1), and a; = 1 without loss 
of generality. When the oscillator reaches (p = 1 aX some moment f, the phase is reset to 
zero and the oscillator produces a pulse signal. If this signal is sent into a delayed feedback 
loop [Fig. [^a)] the emitted pulses affect the oscillator after a delay r at the time instant 
t* = t + T. When the pulse is received, the phase of the oscillator undergoes an instantaneous 
shift by an amount A(y9 = Z{(p{t* — 0)), where Z{ip) is the phase resetting curve (PRC). 
Thus, the dynamics of the oscillator can be described by the following equation [SlinHia]: 

^ = i + z{ip)j2^it-tj-^), ( 1 ) 

b 

where tj are the instants when the pulses are emitted. Note that we adopt the convention 
that positive values of the PRC lead to shorter ISIs. For numerical illustrations we use 
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FIG. 1. (a) Scheme of the model, (b) Shape of the PRC for (7 = 5 (dashed green) and 

q = 28 (solid blue). The stars indicate points with slope —1. (c) Construction of the map Q: 
the new ISI Tj+i depends on the pulse emitted at t = tj-p. (d) Periods T of the RS versus the 
delay r for = 5, obtained from Eq. Q. Solid lines indicate stable RS regimes, dashed unstable. 
Diamonds indicate saddle-node bifurcations. 

Zex{^) ■= 0.1 sin'^ i '^^)) where q controls the steepness of Zex (</?) [see Fig. [^b)]. However, 
onr analysis is valid for an arbitrary amplitnde or shape of the PRC. 

In [H] it was proven that a system with pnlsatile delayed coupling can be reduced to a 
hnite-dimensional map under quite general conditions. To construct the map for system ([^ 
let us calculate the ISI T^+i := tj+i — tj. It is easy to see that T^+i = 1 — Zippj), where 
= ip{t* — tS) = t* — tj is the phase at the moment of the pulse arrival t* = tj_p + r 
|Fig. 0 c)]. Here, P is the number of ISIs between the emission time and the arrival time. 
Substituting tj = tj^p + Tj_p+i + ... + Tj, we obtain the ISI map 

t r*)- 

k=j—P-\-l / 


Tj_|_i — 1 — Z 


T — 


( 2 ) 
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The most basic regime possible in this system is the regular spiking (RS) when the 
oscillator emits pulses periodically with Tj = T for all j. Such a regime corresponds to a 
fixed point of the map ([^ and therefore all possible periods T are given as solutions to 

T = 1-Z{t-PT), (3) 

where P = [t/T], and hence r — PT = r(modT). Figure [^d) shows the period T as a 
function of r for and q = 5. 

To analyze the stability of the RS regime, we introduce small perturbations 6j such that 
Tj = T + Sj, and study whether they are damped or amplified with time. The linearization 
of ([^ in 6j is straightforward and leads to the characteristic equation 

— aX^~^ — aX^~^ — ... — aX — a = 0, (4) 

where a := Z'i^ifj) is the slope of the PRC at the phase = t mod T (cf. [SI [T5]). 

There are two possibilities for the multipliers A to become critical, i.e. |A| = 1. The 
first scenario takes place at a = 1/P when the multiplier A = 1 appears, which indicates a 
saddle-node bifurcation [diamonds in Fig. [^d)]. In general these folds of the RS-branch lead 
to the appearance of multistability and hysteresis between different RS regimes USHH]. 

The second scenario is much more remarkable and takes place at a = —1, where P 
critical multipliers A^ = e*27rfc/(p+i)^ 1 < fc < P, appear simultaneously. This feature is 
quite unusual since in general bifurcations one would not expect more than one real or two 
complex-conjugate Floquet multipliers become critical at once [9]. In the following we study 
this surprising bifurcation in detail and explain why we call it ’’multi-jitter”. 

In order to observe the multi-jitter bifurcation, the PRC Z{ip) must possess points with 
sufficiently steep negative slope Z'{(p) = —1. For instance, in the case Z{(p) = Ze^{(p), 
such points exist for q > q* ^ 27. For such q, two points € (0,1) exist where 

b) = —1 [see stars in Fig. [^b)]. This means that for appropriate values of the delay 
time r, such that r (modT) G a, 4^b}■, it holds a = —1, and the multi-jitter bifurcation 
takes place. Using Eq. (|^ one may determine the corresponding values of r for each possible 
P > 1: 

^A,B = PO- - Z{4a,b)) + 4a,b- ( 5 ) 

Figure shows the numerically obtained bifurcation diagram for q = 28. All values of its 
ISIs Tj observed after a transient are plotted by solid lines versus the delay r. Black lines 
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FIG. 2. (a) ISIs Tj versus the delay r for Q with Z = Zgx and q = 28. Solid lines correspond 
to stable, dashed and dotted to unstable solutions. Black color indicates RS regimes, red stands 
for bipartite solutions (with two different ISIs). Diamonds indicate saddle-node bifurcations and 
stars multi-jitter bifurcations, (b), (c), (d), and (e) are zooms of (a) for P = 1,2,3,4. The lower 
panels show examples of stable bipartite solutions, for which the corresponding sequences of ISIs 
are plotted versus time. Solutions within the same dashed frame coexist at a common value of r, 
which is also indicated by a vertical blue line in the corresponding zoom. 

correspond to RS regimes, while irregular regimes with distinct ISIs are indicated in red 
color. Black dashed lines correspond to unstable RS solutions obtained from Eq. (|^. For 
the intervals r G the RS regime destabilizes and several stable irregularly spiking 

regimes appear. 

Let us study in more detail the bifurcation points t = for different values of P. For 
F = 1, only one multiplier A = — 1 becomes critical. Note that in this case the map ([^ 
is one-dimensional and the corresponding bifurcation is just a supercritical period doubling 
giving birth to a stable period-2 solution existing in the interval r G [Fig. [^b)]. 

For this solution the ISIs Tj form a periodic sequence (0i, © 2 ), where the periodicity of the 
sequence is indicated by an overline. It satisfies 

Q 2 = 1 — Z {t — ©i) and ©1 = 1 — Z (r — © 2 ). (6) 

For P > 2, P multipliers become critical simultaneously at r = and the RS solution 
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is unstable for r G Numerical study shows that various irregular spiking regimes 

appear in this interval. We observe solutions, which have ISI sequences of period P + 1 but 
exhibit only two different ISIs in varying order [see Fig. bottom]. As a result, each solution 
corresponds to only two, and not P + 1, points in Figs. |^a),(c)-(e). In the following we 
call such solutions “bipartite”. For larger P, a variety of different bipartite solutions with 
(P + l)-periodic ISI sequences can be observed in r G The stability regions of 

these solutions alternate and may overlap leading to multistable regimes [see Appendix, 


Figs, A.1 - A.4 


The bipartite structure of the observed solutions can be explained by their peculiar 
combinatorial origin. Indeed, all bipartite solutions can be constructed from the period- 


2 solution(0i, © 2 ) existing for P = 1. Consider an arbitrary (P + l)-periodic sequence of 


ISIs (Ti,T 2 , ...,Tp+i), where each Tj equals one of the solutions © 1^2 of (|^ for some delay 
r = To G ['T’^jT'^j- Let ni > 1 and n 2 > 1 be the number of ISIs equal to ©1 and to ©2 
respectively. Then it is readily checked that the constructed sequence is a solution of ([^ at 
the feedback delay time 


Pil,n2 — T'o + (^1 ~ 1) ©1 + (^2 ~ 1) 02- 


(7) 


Red dotted lines in Figs, [^c), (d), and (e) show the branches of bipartite solutions con¬ 
structed from ([^ with P = ni -t- n 2 — 1 = 2, 3,4. Note that these solutions he exactly on 
the numerical branches which validates the above reasoning. However, some parts of the 
branches are unstable and not observable. Since each bipartite solution corresponds to a 
pair of points Ti) and {Tn^,n 2 ,T 2 ), solutions with identical ni and n 2 correspond to the 

same points in the bifurcation diagrams in Fig. For instance, when P = 3 the branches 
corresponding to the solutions (© 1 , © 2 , © 1 , © 2 ) = (© 1 , © 2 ) and (© 1 , © 1 , © 2 , © 2 ) he on top of 
each other. 

Let us estimate the number of different bipartite solutions for a given P G N. The number 
of different binary sequences of the length P -|- 1 equals 2^^^. Subtracting the two trivial 
sequences corresponding to the RS one gets 2^+^ — 2. Disregarding the possible duplicates 
by periodic shifts (maximally P -|- 1 per sequence) one obtains an estimate for the total 
number Mp of bipartite solutions for a given value of P as 


A/p P 


2 ^+^ - 2 

Xp+iy 


( 8 ) 
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Notice that all these bipartite solutions exist at the same value of P but for different ranges 
of the delay r. Nevertheless, all emerge from the RS solution in the bifurcation points 
To see this, let us consider the limit tq —)■ t\. In this case the ISIs 0i(ro) and 02 (t'o) tend 
to the same limit Ta = 1 — Z {iPa) which is the period of the RS at the bifurcation point. 
Then, Q converges to + (ni + 77-2 — 2)Ta = Ta, while all bipartite solutions converge to 
the RS with period Ta- 

Thus, all bipartite solutions branch off the RS in the bifurcation points This hnding 
is clearly recognizable in Figs. |^c)-(e), where stars indicate the multi-jitter bifurcation 
points. Numerical simulations show that many of the bipartite solutions stabilize leading to 
high multistability. In particular, we observe that all bipartite solutions with same values 
of rii and n 2 exhibit identical stability. This emergence of numerous irregular spiking, or 
jittering, regimes motivates the choice of the name “multi-jitter bifurcation”. 

High multistability is a well-known property of systems with time delays. A common 
reason is the so-called reappearance of periodic solutions [18]. This mechanism may cause 
multistability of coexisting periodic solutions, whose number is linearly proportional to the 
delay. Due to multi-jitter bifurcation, multistability can develop much faster, since the 
number of coexisting solutions grows exponentially with the delay [cf. Eq. ([^]. This suggests 
that the underlying mechanism is quite different. 

Irregular spiking regimes similar to the ones described here were reported previously for 
systems exhibiting a dynamical “memory effect”, where the effect of each incoming pulse 
lasts for several periods [T9|. In system ([^, however, the effect of a pulse decays completely 
within one period. Therefore the origin of jittering must be different. In fact, it relies on 
another kind of memory, which is provided by the delay line and stores the last P ISIs. This 
memory preserves fundamental properties of time, which are responsible for the degeneracy 
of the multi-jitter bifurcation. To explain this let us consider (|^ as a P-dimensional mapping 
(Tj-P+i, ...,Tj) I— )■ (Tj_p+ 2 , ...,Tj+i). Disregarding the calculation of the new ISI T^+i as in 
(|^, all the map does is to move the timeframe by shifting all ISIs one place ahead. The rigid 
nature of time allows no physically meaningful modihcation of this part of the map which 
could unfold the degenerate bifurcation. Moreover, the new ISI T^+i depends exclusively 
on the sum of the previous ISIs which has the effect that all past intervals have an equal 
influence on the new ISI regardless of their order. As a consequence the combinatorial 
accumulation of coexisting solutions with differently ordered ISIs is generated. 


Besides the delayed feedback, another essential ingredient for the mnlti-jitter bifnrcation 
is the existence of points where the PRC fnlhlls Z'{ip) < — 1. Since the PRC is a characteristic 
that can be measnred for an arbitrary oscillator [S] , this condition gives a practical criterion 
for the occnrence of jittering regimes. In this context it is worth noticing that the condition 
Z'{ip) < — 1 is eqnivalent to the non-monotonicity of the system’s response to an external 
pulse, i.e. there are such phases (pi < ip 2 that a pulse can reverse their order as (pi + Z{(pi) > 
Lp 2 + Z{lp 2 ). Note that for a smooth one-dimensional system as Q the reversal of phases is 
only possible if the feedback takes the form of 5-pulses. With pulses of hnite duration two 
continuous orbits connecting the different phase values before and after the pulse cannot 
cross each other. This prevents a reversal of phases. However, for oscillators with a phase 
space of dimension larger than one the phase points (fi and ip 2 can exchange their order 
without necessitating the orbits to intersect. 

In order to evaluate the practical relevance of the theory described above we consider two 
realistic systems: (i) an electronic implementation of the FitzHugh-Nagumo oscillator |20l- 
l23] with time-delayed pulsatile feedback, and (ii) a numerically simulated Hodgkin-Huxley 
model | 23 ] with a delayed, inhibitory, chemical synapse projecting onto itself HSIIIT]. A 
detailed description of the systems is given in the Appendix [Secs. and . In both cases 
the measured PRC exhibits parts with slope less than —1 [ see Figs. B. I|and [CT] . Therefore 
the existence of multi-stable jittering can be conjectured on the basis of our results for 
system Q. Figure presents experimentally (for the FitzHugh-Nagumo oscillator) and 
numerically (for the Hodgkin-Huxley model) obtained bifurcation diagrams showing ISIs for 
varying delays. Both systems clearly show that a stable RS solution destabilizes closely to 
the multi-jitter bifurcation points. Where the RS regime is unstable, the system switches 
to irregular spiking, and we mainly observe {P + l)-periodic bipartite solutions. The insets 
in the lower part of Fig. |^a) show two such period-5 bipartite regimes of the FitzHugh- 
Nagumo oscillator. Note that both of them coexist at r = 126 ms, which illustrates the 
multistability of the system. Similarly, two period-4 bipartite regimes coexisting for r = 60 
ms for the Hodgkin-Huxley model are shown in|^b). 

Bipartite solutions are a basic form of jittering both in phase-reduced models and re¬ 
alistic systems. However, beyond the multi-jitter bifurcation the bipartite solutions may 
undergo subsequent bifurcations. In system ([^, we observed higher-periodic, quasiperiodic 
and chaotic regimes for larger steepnesses of the PRC {q > 70). Similar regimes were also 
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FIG. 3. Bifurcation diagrams for (a) an electrically implemented FitzHugh-Nagumo system and 
(b) a simulated Hodgkin-Huxley neuron model. RS solutions are indicated by black dots, solutions 
with more than one ISI by red dots. The branches of RS solutions obtained from the measured 
PRCs and Q are plotted as dashed grey lines. Multi-jitter bifurcations are indicated by star¬ 
shaped markers. The values of the delay r are given in ms. The panels below the diagrams show 
examples of stable bipartite solutions (ISIs versus time). The values of r for which the solutions 
exist are also indicated by vertical blue lines in the corresponding diagrams. 

found for the Hodgkin-Huxley model. An example showing aperiodic jittering is shown in 
Fig. [^b) for r = 59.5 ms. In the Appendix we show more examples of aperiodic jittering 
[see Fig. |D.1| . 

To conclude, in a phase oscillator with delayed pulsatile feedback Q we discovered a 
surprising bifurcation leading to the emergence of a large number (~ expr) of jittering 
solutions. We showed that this multi-jitter bifurcation does not only appear in phase- 
reduced models, but also in realistic neuron models and even in physically implemented 
electronic systems. These hndings support our theoretical results and provide motivation 
for a deeper study of the multi-jitter phenomenon. 

The possibility of jittering depends on the steepness of the PRC which is an easily mea¬ 
surable quantity for most oscillatory systems [5]. Thus, our hndings provide an easy criterion 
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to check for the existence of jittering in a given system. This may prove useful in a vari¬ 
ety of research areas, where pulsatile feedback or interactions of oscillating elements takes 
place. For instance, this might be one of the mechanisms behind the appearance of irregular 
spiking in neuronal models with delayed feedback [25] and timing jitter in semiconductor 
laser systems with delayed feedback [26] . For applications which exploit complex transient 
behavior such as liquid state machines IZ71 the high dimension of the unstable manifold at 
the bifurcation can be interesting. Furthermore, in view of the possibility of a huge number 
of coexisting attracting orbits beyond the bifurcation the system can serve as a memory 
device by associating inputs with the attractors to which they make the system converge. 
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Appendix A: Maps and basins for P = 1,2,3,4 


In this section we illustrate the ISI maps for cases P = 1,2, 3,4. In all cases we use 
^ex(p) := 0.1 sin'^ with q = 28. 

P = 1. In this case the map is one-dimensional: 


Tj+i = 1 -Z(t-Tj). 


(A.l) 


The dynamics of the map can be illustrated on a coweb diagram which is depicted in 
Fig. 1^ for T = 1.5. For this value of the delay, the only attractor is a stable period 2 


solution [cf. Fig. 2 of the main text]. 



WM 


FIG. A.l. Cob-web diagram of (A.l) for r 


1.5. The black solid dots correspond to a stable 


period 2 solution and the white hallow dot corresponds to unstable regular spiking. 
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P = 2. In this case the map is two-dimensional and reads 


— 1 — Z{t — Tj — Tj_i). 
This can be equally rewritten in the vector form: 


(A.2) 


/ Ti 1 

u( 

1 T2 i 



T2 

1 -Z(r-T2 -Ti) 


Figure A.2 shows attractors and their attraction basins of the map (A.2) for r = 2.414. 


For this value of the delay we observe a coexistence of a stable regular spiking solution and 
a stable irregular period-3 solution of the form (0i, ©i, © 2 ) [cf. Eq. [^of the main text]. 

0.98| 



0.96 


0.94 


0.9 0.92 0.94 0.96 0.98 


FIG. A.2. Phase plane of (A.2) for r = 2.414. The black solid dots correspond to a stable period 


3 solution and the white hallow dot corresponds to stable regular spiking. The blue region is the 
attraction basin of the period-3 solutions, red of regular spiking. 
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P = 3. In this case the map is three-dimensional and reads 


Tj+i — 1 - Z {t - Tj - Tj_i - Tj_2), 


(A.3) 


or, written in an equivalent vector form 

/ 


f T^\ 

T2 

\Tsj 


-)■ 


T2 

Ts 




^ 1 — Z (r — T 3 — T 2 — Ti) j 


For T = 3.37 three different stable bipartite solutions of this map coexist: (0i, ©i, ©i, © 2 ), 


(© 1 , © 1 , © 2 , © 2 ), and (© 1 , © 2 ). Figure A.3 depicts the basins of attractions conhned to the 
two-dimensional plane 


P = {(Ti,T 2,T3)| Ti,T2e [0.9,0.98],T3 = Ti} (A.4) 

of the three-dimensional phase space. 



0.9 0.92 0.94 0.96 0.98 


T, 


FIG. A.3. Planar section of the attractor basins of the map (A.3) by the plane (A.4) for r = 3.37. 


The corresponding attracting bipartite solutions are depicted in the right part. 
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P = 4. In this case the map is four-dimensional and has the form 


Tj+i — 1 — Z (t — Tj — Tj_i — Tj_2 — Tj_^). (A.5) 

Here we omit the vector form for brevity. For r = 4.32 four different stable bipar¬ 
tite solutions coexist: (01, 01, 02, ©2, ©2), (01,02,01,02,02), (01,01,01,02,02), and 
(01, 01, 02, 01, 02)- An intersection of the attractors basins with the 2 -dimensional plane 


H 


{(Ti,T2,T3,T4)|Ti,T 2 e [0.9,0.98],T3 = Ti,T4 = T2} 


is shown in 


Fig. 





(A.6) 


FIG. A.4. Planar section of the attractor basins by the plane (A.6) of the map (A.5) for and 


r = 4.32. The corresponding attracting bipartite solutions are depicted in the right part. 
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Appendix B: Electronic FitzHngh-Nagnmo oscillator 


The circuitry of the electronic FitzHugh-Nagumo oscillator as used in the experiment 
[cf. Fig. I^a) of the main text] is depicted in Fig. B.l[ a), see Refs. [221 123] for details. 
Here, R = Ikhl, C = 47nF, L = 103.4H, Pin is an input from the delay line, and F{u) = 
au{u — uq){u + uo) is the cur rent-volt age characteristic of the nonlinear resistor with a = 
2.02 X 10 ^ and uq = 0.82V. The autonomous oscillations have period T 30ms. 

The delay line is realized as a chain of monostable multivibrators. A pulse of the amplitude 
A = 5V and duration 0 = 0.42ms is delivered with a delay r each time the voltage u reaches 
the threshold Uth = —0.7V. 

The noise level in the circuit was ~ —40dB (this means fluctuations of ~ 20mV for a 
signal amplitude of ~ 2V). 

The measured phase resetting curve for the given parameters is depicted in Fig. [B.lK b). 
The interval where the PRC slope is less than minus one is highlighted in red. 




FIG. B.l. (a) Circuitry of the Fitzhugh-Nagumo oscillator, (b) The measured PRC of the 
FitzHugh-Nagumo oscillator. The interval with the slope less than minus one is highlighted in 
red. 
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Appendix C: Hodgkin-Huxley model 

The periodically spiking Hodgkin-Huxley neuron model, which was used for the numeri- 
calresults in Fig. [^b) of the main text, is given by the following set of equations [muziEi] 


CV{t) = 1- gMam^h{V{t) - VNa) - 9Kn{V{t) - Vk) 

- gi{V{t) - Vi) - K{V{t) - Vr)s{t - r), (C.l) 

m{t) = am{V{t)){l -m{t)) - 

h{t) = ah{V{t)){l - h{t)) - (]h{V{t))h{t), 
n{t) = an{V{t)){l -n{t)) -/3„(H(t))n(t), 
s{t) = 5(1 - s(t))/(l + exp(-l/(t))) - s(t). 


where V{t) models the membrane potential, amiV) = (O.IH -|- 4)/(l — exp(—O.IH — 4)), 
l3m{V) = 4exp((-l/-65)/18), ah = 0.07exp((-l/- 65)/20), l3h{V) = 1/(1exp(-0.1H- 
3.5)), an{V) = (O.OIH + 0.55)/(l - exp(-0.1H - 5.5)), /?„(H) = 0.125exp((-H - 65)/80), 
C = 1/iF/cm^, / = lOpA/cm^, g^^a = 120mS/cm^, V^a = 50mV, gK = 36mS/cm^, Vk = 
—77mV, gi = 0.3mS/cm^, Vi = —54.5mV, K = —65mV, and n = 0.38mS/cm^. 

Figure [CT| shows the PRC of (C.l), which was measured by replacing the delayed feedback 
by an external stimulation line through which singular synaptic pulses sampled from the 
same system were applied at different phases. 



FIG. C.l. PRC of the Hodgkin-Huxley model (C.l). The region, where Z' (<p) < —1 is indicated 
by red. 
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Appendix D: Emergence of chaotic jittering in the Hodgkin-Huxley model 


Figure D.l illustrates the emergence of chaotic jittering states for increasing feedback 
strength k. It shows three different trajectories for r = 59.3 and k = 0.38, 0.4, 0.41. In plot 
(a), the ISIs of a quasiperiodic solution are shown for c = 0.38. A black dot is placed at 
(tj,Tj), where tj is the moment when the j-th ISI ends and Tj is its duration. Subsequent 
dots are joined by a blue line. The sequence is contained in a torus in the phase-space, whose 
projection to the (T^,T^+ij-plane is shown in plot (b). For each Tj from the sequence of 
approximately thousand observed ISIs, a blue dot was placed at the coordinates {Tj,Tj^i). 
Plots (c) and (d) depict a solution close to the onset of chaotic jittering for c = 0.4. The 
corresponding Lyapunov exponent is positive but small [see Fig. D.2[ b)]. A more pronounced 
chaos is exhibited by the solution existing at c = 0.41, which is shown in plots (e) and (f). 
Note that the emergence of chaos is accompanied by a loss of smoothness of the torus 
consisting of the quasiperiodic trajectories [28] . 

Figure [D^ shows the numerically calculated Lyapunov exponents (LE) for different values 
of the feedback strength k. For each value of k the four largest exponents are shown. 
Starting from r = 58 on the RS solution with maximal period [cf. Fig. [^b) for k, = 0.38] 
each computation for r G [58, 60.5] was initialized with a solution on the previous attractor 
as initial data. Note that there exists always one LE which has real part zero, since the 
corresponding attractors are not steady states. For each value of r, the point where the 
multi-jitter bifurcation occurs can clearly be recognized as all four depicted LE approach 
zero nearly at the same value of t = ta {ta ~ 58.75 for k = 0.38, ta ~ 58.4 for n = 0.4, 
and Ta ~ 58.3 for k = 0.41). Dynamics found in the depicted range beyond these points, 
i.e. for Ta < T < 60.5, are irregularly spiking regimes. For k = 0.38, two LE are zero in the 
interval r G [59.25,59.8], which indicates a torus, i.e. quasiperiodic behavior as illustrated 
in Fig. D.l[ a). For all other values of r in the depicted range, we observe periodic solutions 
which exhibit ISI sequences of period 4 [see Figj^b) of the main text]. When the feedback 
strength is increased to k = 0.4 the dynamics in the corresponding interval r G [59.2, 59.5] 
become weakly chaotic. For k = 0.41 the LE become larger in the interval r G [58.75, 59.2] 
which indicates a more pronounced chaos. Note that even for feedback strengths where 
chaotic jittering is observed, periodic solutions still exist at other values of r. 
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FIG. D.l. Emergence of chaotic jittering in (C.l). The plots (a), (c), and (e) show ISI sequences of 


trajectories observed at the indicated feedback strengths k = 0.38, 0.4, 0.41. (a),(b): Quasiperiodic 
jittering for c = 0.38; (c),(d): weakly chaotic jittering for c = 0.4; (e),(f): clearly recognizable 
chaotic jittering for c = 0.41. 
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FIG. D.2. Four largest Lyapunov exponents of (C.l) with r G [58,60.5] and feedback strengths 
(a) K = 0.38; (b) k = 0.4; (c) k = 0.41. Arrows indicate the delays, for which the trajectories in 


Fig. D.l are shown. 























